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O ■ Abstract 

■ We show faster algorithms for solving regression problems based on estimating statistical leverage 
(N. scores. A growing number of applications involve large, sparse, n x d matrices A where n ^ d. For 

many of these applications the more expensive operations involve the d x d matrix A"^A. When n is 
much larger than d, the running time bottleneck is in the cost of computing the matrix A^A, rather 

■ than the more expensive operations on this smaller matrix. 

, Recent works by Clarkson, Drineas, Magdon-Ismail, Mahoney and Woodruff led to algorithms that 

approximate A'^A in 0(n(i log n + poly((i)) and 0(nnz(A) + poly((i)) time respectively. When the size 
, of A is much larger than the cost of more expensive operations involving the smaller dx d matrix, these 

algorithms offer significant savings in running time. 

We give alternate approaches for approximating A^A based on techniques originally developed for 



^ ■ spectral sparsification. Our algorithm finds a matrix B consisting of a small number of scaled rows of A 

CO , such that with high probability ||Ax||2 = (1 ± e)||Bx||2 for all vectors x e The running time of our 

algorithm can be bounded by 0{nnz{A) + (i"+"e~^) for any constant a > 0. The key to our approach 
, is to find a sequence of estimates of A^A and gradually improve the approximations. 

1 Introduction 

. Least squares calculations and singular vector /value computations are among the most common compu- 
[ tational operations involving matrices. These routines are used as core steps in clustering, learning, and 
K*" ■ more recently algorithmic graph theory. For many operations involving the matrix A, the computationally 
^ . more expensive steps are often performed on its product, A A. One of the simplest settings where this 
5h [ occurs is with over-determined systems of linear equations of the form Ax = b. Multiplying both sides 
- - ■ by A gives A Ax = A b, and when A has full column rank, multiplying both sides by (A A)~ gives 
X = (A^A)-iA^b. This can be shown that this choice of x minimizes 1 1 Ax — b| I2, a common objective in 
least squares regression [Str93]. 

Many of the operations involving A"^A, such as inverse or eigenvector computation, take 0{d^) time 
where w is the matrix-multiplication constant. When the number of rows, n, is comparable to the number 
of columns, d, A'^A can be computed at a fraction of the cost of these more expensive operations. Also, 
for iterative methods, the linear operator corresponding to A^A can be implicitly applied in 0(nnz(A)) 
time via. two matrix-vector multiplies involving A and A"^, where nnz(A) is the number of non-zero 
entries. As a result, having access to the linear operator corresponding to A"^A is taken for granted in 
most algorithms. However, many applications involve matrices that are tall and skinny, where n is much 
larger than d. Here, the d"^ entries in A"^A is smaller than the original number of non-zero entries in A 
itself. Once n is more than d^, the cost of computing A"^A is the main running time bottleneck both in 
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theory and in practice. Direct evaluation of A A takes 0(nnz(A)(i) time when A is sparse, or 0{nd^~^) 
time when A is dense. These instances are sufficiently common that practical speedups for finding QR 
factorizations of tall and skinny matrices have been studied in the distributed [SLHDIO, ACD+10] and 
MapReduce settings [CGll]. 

A more direct way of obtaining speedups is to find faster algorithms that give approximations to A"^ A. 
The error of approximation between B^B and A"^A can be measured either additively or multiplicatively. 
A multiplicative guarantee is more appealing for several reasons. It can be viewed as relating |[Axj|2 to 
||Bx|j2 for all vectors x. This means that B can be directly used in place of A for the overconstrained 
system of linear systems mentioned earlier. Furthermore, this guarantee is same as the requirement of 
usmg B^B as a preconditioner for A^A in iterative methods (see e.g. [Axc94, Saa96]). This means B^B 
can be used in place of A^A to obtain higher accuracy solutions with only a constant factor overhead. 
This connection is crucial in practice, and has been used to give speedups on overconstrained systems over 
the solver in LAPACK [AMTIO]. 

One of the first studied approaches for approximating A"^ A is to find a matrix B consisting of a subset 
of rescaled rows of A, and then compute B^B. If B has fewer rows, B^B can be evaluated faster, and it's 
also possible for it to remain close to A"^A. To our knowledge it was first shown by Drineas, Mahoney, 
and Muthukrishnan that this type of algorithms can be used to speed up L2 regression [DMM06], albeit 
in restricted settings. For general matrices, Magdon-Ismail first showed that B can be found in time 
faster than directly computing the product A^A [MHO]. Drineas et al. [DMIMW12] gave a much faster 
algorithm for row sampling that runs in 0{ndlogn + poly((i, e)) time. Recently, Clarkson and Woodruff 
[CW12] gave algorithms for computing least squares regression and low rank approximations in input 
sparsity, or 0(nnz(A) + poly((i, e)) time. 

These algorithms have in common that they apply a sequences of transforms to A oblivious to its 
structure, and arrive directly at the high-accuracy (l±e) approximations. For non-uniform A concentrated 
in a few rows, it's clear that these rows are more important. In fact, non-uniform sampling was one of 
the earliest studied methods [DMM06], and remains the method of choice if key sampling probabilities 
are provided [DMIMW12]. This paper gives an adaptive algorithm based on row sampling that runs in 
0(nnz(A) + time. Our techniques are motivated by the connection between row sampling and 

computing statistical leverage scores, and are motivated by works by Spielman and Srivastava on spectral 
sparsification of graphs [SS08]. In this setting, A is the edge- vertex incidence matrix of a graph, and 
subsequent extensions of their work showed that sampling probabilities can be bounded by stretches w.r.t. 
a spanning tree [KLP12]. Our approaches generalize this notion to stretch w.r.t general matrices, and we 
denote these matrices that we use to measure stretch as bases. These generalized stretches can be used as 
sampling probabilities, and we will show that similar matrices can be used in place of each other as bases. 
This allows us to compute approximate stretches and row samplings in an iterative manner, giving the 
following result: 

Theorem 1.1 Given a n x d matrix A along with failure probability 6 = d^^ and allowed error e. For 
any constant a > 0, we can find in 0{nnz{A) + d'^'^'^ e~~'^) time, with probability at least 1 — 6, a matrix 
B consisting of 0{{d^^" + ""J^a^^ )^~'^) rescaled rows of A such that I|Aa;||2 = (1 =be)|jBa;||2 for all vectors 

The organization of the paper is as follows. Section 2 describes notations used throughout the paper, 
reviews known results that describe key sampling probabilities, and defines generalized stretch. These 
results, which are mostly present in earlier works, give that crude estimates of A"^A can lead to better 
ones. Then we show two iterated algorithms for row sampling in Sections 3 and 4. Section 3 gives an 
algorithm that iterates on projected versions of A, giving the result as stated in Theorem 1.1. Section 4 
describes an alternate algorithm that iterates on an artificially perturbed version of A, giving an algorithm 
with logarithmic dependency the condition number of A^A. 
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2 Background 

In this section we outline key results from previous works on row sampling. These works show that statis- 
tical leverage scores are closely associated the probabilities needed for row sampling, and give algorithms 
that efficiently approximate these values. In addition, we formalize the observation implicit in earlier works 
that both the sampling and estimation algorithms are very robust. The high error-tolerance of these algo- 
rithms makes them ideal as core routines to build iterative algorithms upon. Most of the results mentioned 
in this section are direct adaptations from earlier works, where they're presented in much more detail. For 
completeness, we give their proofs in Appendix A. 

Most of the presentation will use standard linear algebraic notations. We use ||x||p to denote the Lp 
norm of a vector. The two values of p that we'll use are p = 1 and p = 2, which correspond to ||x||i = \xi\ 

and ||x||2 = Y^^j xf. For two vectors x and y, x > y means x is entry- wise greater or equal to y, aka. 

Xi > Yi for all i. Given a matrix A, we use Aj^;, or to denote the i^^ row of A. Note that a, is a 
row vector of length d. We will also use its Forbenius norm, HAHj? to denote the 2-norm of all its entries. 

Formally, ||A||| = ^^^Af". 

The spectral theorem states that any symmetric matrix C has a full complement of orthonormal 
eigenvectors. Let them and their eigenvalues be (ui, Ai), (u2, A2), • • • (u„, A„). Here Uj is a length d column 
vector, Xj is a scalar and they satisfy A"^Auj = Xjuj. Without loss of generality we assume that AjS are 
in increasing order: Ai < A2 < • • • < A^. Then if A^ is the first non-zero eigenvalue, we have: 

n 
j=k 

We will use Amax to denote the largest eigenvalue, and Amin to denote the smallest non-zero eigenvalue. 
Suppose Afc is the first non-zero eigenvalue, the pseudoinverse of C, can be defined as: 

j=k ^ 

This is a linear operator that preserves the null space of C, while acting as its inverse on the rank space. A 
matrix C is positive semi-definite if all its eigenvalues are non-negative, or equivalently x-^Cx > for all 
vectors x. Similarities between matrices is defined using a partial order on matrices. Given two matrices 
A and B, A ^ B denotes that B — A is positive semi-definite. 

We will use K to specifically to denote A^A. An alternate view of row sampling is that we're given n 
PSD matrices Ki . . . K„ where Kj = a^a^; and the goal is to find a small subset of them such that their 
weighted sum is close to K. Furthermore, we can assign weights of zero to elements that are not selected, 
giving a weight vector w G 3^". The formal requirement of row sampling can be described as: 

• For some parameter k we have K ^ E"=i WiK.i ^ kK 

• \{i : Wi ^ 0}\, the support of w, is small. 

The work by Glarkson and Woodruff [CW12] obtained approximations to A by applying a (randomized) 
linear transformation S to get B = SA. In our case, such a matrix can be constructed by letting each row 
correspond to a non-zero weight in w. This row in S then has one non-zero entry equaling to the square 
root of the weight. This is similar to the method used by Glarkson and Woodruff in that each column has 
at most one non-zero entry, but has the additional property that each row has at most one non-zero entry, 
that's also non- negative. 
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2.1 Row Sampling Using Statistical Leverage Scores 

Matrix concentration bounds can be viewed as generalizations of single-variable concentration bounds. To 
our knowledge, the first nearly-tight bounds that encapsulates row sampling was given by Ahlswede and 
Winter [AW02]. Similar bounds were given as part of a much larger set of results by Rudelson and Ver- 
shynin [RV07], although simplifications of their proofs are similar to previous works [Vcr09, Harll]. These 
results identified statistical leverage scores as the crucial numbers corresponding to sampling probabilities. 
Formally, the statistical leverage score of a row of A, a^ is defined as: 

Ti = ajK+af 

Where is the pseudoinverse of K. 

Although these values have been studied in statistics, their use in algorithms is more recent. Spielman 
and Srivastava used these bounds in their spectral sparsification algorithm [SS08]. However, their algorithm 
was restricted to the setting where A is the edge-vertex incidence matrix of some graph. Works by Drineas, 
Mahoney, Avron and Toledo generalized this sampling scheme to arbitrary PSD matrices [DM10, ATll]. 
For our algorithms, we have crude approximations of leverage scores, but only need to bound the number 
of samples. Such extensions are immediate from the presentations given in [Ver09, Harll]. We use the 
statement obtained by combining Theorems 5.4. and 6.2. from [ATll], as it gives guarantees for sampling 
with upper bounds of leverage scores. 

Lemma 2.1 (Theorem 5.4, 6.2 from [ATll]) Let A be a n x d matrix and K be its inner product, K = 
a'' A. Suppose f is a length n vector such that fj > Tj. Associate the probability with every row: 

m=u^ (2.1) 

Irlli 

and let Ti . . . T^j be i. i. d. random matrices defined by: 

T,=pj^Kj^ (2.2) 

Where Ji . . . Jm are random integers between 1 and n which take value i with probability pi. Then there 
exist an absolute constant C such that given any error parameter e, choosing M = C||f||i log(||f||i)/e^ 
gives: 



^ (l + e)K >1 



poly{M) 

Furthermore, this process, which we denote as Sample(A, f, e) can run in 0{n + Mlogn) time. 

Note that the total number of sample depends on the sum of the upper bounds. In the case we use the 
exact scores, the following fact implies about O(dlogd) samples suffice. 

Fact 2.2 (see e.g. [SS08j) 

n 

y^^Tj = r < d 
1=1 

The sampling process described in Lemma 2.1 is with replacement, and is similar to the framework 
described in [AW02, RV07, Harll]. The running time can be obtained by computing partial sums of f and 
binary search for the entry where the cumulative probability meets our threshold. Ipsen and Wentworth 
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showed that a variety of different samphng methods lead to similar results both theoretically and experi- 
mentally [IW12]. Since fj only needs to be upper bounds for Tj, we can also use poor approximations of 
leverage scores by scaling them up by the approximation factor. This scaling leads to an increases in the 
number of rows sampled by the same factor as the approximation ratio. However, it allows us to obtain a 
much more accurate estimate of A compared to the accuracy of the approximate statistical leverage scores. 

2.2 Generalized Stretch and Their Properties 

As our algorithm is iterative and its output is approximate, we will not be able to compute K exactly. 
Therefore, all the statistical leverage score estimations in our algorithms will be done using an approxima- 
tion K in place of K. As a result, we will incorporate the matrix that we compute these scores by, K as 
part of our notation. This is a common notion used in the design of combinatorial preconditioners, where 
K is often chosen to be a simpler matrix than K. In the tree setting, Spielman and Woo [SW09] showed 
that when A is an edge-vertex incidence matrix and K corresponds to a spanning tree, ajKa?" equals to 
the combinatorial stretch of the edge i w.r.t. the tree. This motivates us to use generalized stretch to 
describe this term for arbitrary matrices. Suppose we're given a factorization of K as B"^B, we will use 
5T7^B(aj) to denote: 

5r7^B(ai) := aiK+a, = a,(B^B)+a, (2.3) 

Note that under this notation, the original definition of statistical leverage score of row i, ti equals to 
STTZAi^i)- We will also use this generalized notion to denote leverage scores from this point and on, but 
still use f to denote the approximate bounds of leverage scores computed by our algorithms. 

Instead of measuring the stretch of a single row of A, we can also incorporate the stretch of a subset 
of the rows into our notation. If A is a matrix with n rows, we use the notation 5T7^b(A) to denote: 

n 

5r7^B(A) = J^5r7^B(ai) (2.4) 

1=1 

The following fact is also crucial for our analysis: 

Fact 2.3 Let K = B. The generalized stretch of the i^^ row w.r.t B equals to its L\ norm under the 
transformation K^^^"^ : 

STli-B\Od) =\\J^ O-i II2 = \ \0-iK II2 

and the total stretch of all rows is 



sttzb{a) =\\it''^A'\\l = WBit^'^F 

— hl/2 

Where K is the 1/2 power of the pseudo-inverse of K and \ \ ■ \ \f denotes the Forhenius norm. 

Computing stretch using a different matrix, 5T7^b(A), is a standard step when analyzing graph based 
preconditioners, and done implicitly in [KLP12]. It was shown by Avron, Maymounkov and Toledo that 
similar matrices can be used in place of each other in least squares regression [AMTIO]. The same also 
holds for generalized stretch, where similar matrices can be used to approximate stretch up to the same 
factor as their difference. This allows us to define the notion of S-stretch base. 

Definition 2.4 A matrix B is a S -stretch base for A if: 

1. STlZBifh) > Ti 
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2. STTZb{A) < S 



This definition is similar to the objective of the column subset selection problem studied by Avron 
and Boutsidis [AB12]. Its two conditions are consequences of the requirements on sampling probabilities 
described in Lemma 2.1. Condition 1 implies that 5T7^A(aj) can be used as sampling probabilities; while 
condition 2 means if k = poly((i), the total number of rows needed for a constant error can be bounded 
by 0{ndlog d). Works on spectral sparsifiers of graphs showed that matrices similar to A, B, can still lead 
to sampling probabilities [SS08, KLP12]. One sufficient condition is that B"^B is spectrally close to A"^A. 
It can be formalized using Definition 2.4 as follows: 

Lemma 2.5 For fixed A and K = A, if K = B satisfies: 



Then B is a K-stretch base for A. 

Using diff'erent bases play crucial roles in both of our algorithms described in Sections 4 and 3. 
2.3 Estimation of Generalized Stretch 

Fact 2.3 shows that gives that STlZxi^ii) is the L2 norm of a vector. Combining this norm using the 
Johnson-Lindenstrauss transform (e.g. [AchOl]) the leads to faster approximations of all leverage scores. 
This was first used by Spielman and Srivastava to approximate effective resistances [SS08]. Direct appli- 
cations of this reduction, or similar algorithms given in [DMIMW12] lead to a running time of at least 
J7(nnz(A) log n). However, the definition of generalized stretch are only one sided bounds pointwise, which 
means our estimates can also have large, one sided error. We can use fewer vectors in the projection, and 
scale up the results to correct potential underestimates. This trades off the coefficient on the leading term 
at the cost of more elements sampled. This way of remedying one-sided error in the Johnson-Lindenstrauss 
transform was first shown in [KLP12] and relies on a version of the projection theorem with fewer vectors. 
We use the following version given as Lemma 7 in [IM98]: 

Lemma 2.6 (Lemma 7 from [IM98]) Let u he a unit vector in . For any given positive integers k, let 
Ui, . . . ,Uk be random vectors chosen independently from the v-dimensional Gaussian distribution N'^(0, 1). 
For Xi = u^Ui, define W = W{u) = (Xi, ...,Xk) and L = L{u) = \\W\\l. Then for any p > 1: 

1. E{L) = k 

2. Pr{L > pk) < 0(A;)exp(-|(/3- (l + ln/3))) 

3. Pr{L < k/P) < 0(A;)exp(-|(/3-i - (1 -ln/3))) 

~ +1/2 

The other observation needed to apply this lemma is that K from Fact 2.3 can be replaced by any 

matrix whose product with its transpose equals to K. As K = B"^B, = K^B-^BK"*" = (BK^)-^BK'^. 

~ +1/2 

Using this matrix in place of K avoids computing the 1/2 power of K. Pseudocode of our estimation 
algorithm is shown in Algorithm 1, while the error analysis is nearly identical to the ones given in Section 
4 of [SS08] and Section 3.2. of [DMIMW12]. 

Lemma 2.7 For each i, with probability at least 1 — 5 we have: 



-K<k< K 



K 
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Algorithm 1 Algorithm for Upper Bounding Stretch 



ApproxStr(A, B, /3, 6) 

Input: Matrix A, approximation B such that B"^B is a K-stretch base for A, value k, parameter /3 indicating 
allowed estimation error. Failure probability 6. 

Output: Upper bounds for stretches of rows of A measured w.r.t. K: fi . . . f„. 

Compute K = B'^B 

Compute K"*" 

Let k = 0(log^(l/(5)) 

Let JJ he a k X n matrix with each entry is picked independently from A^(0, 1) 
Let n be |||UAK^af ||^. 
return f 



T 



Proof We can rewrite K as K A AK to get: 



5T7^J>(aj) =aiK a^ 



=aiK KK a/ 
=aiK^(A^A)K^a; 
=aiK"^A^AK"^af 
= l|AK+aflli 



(2.5) 



By Lemma 2.6 Part 3 we have: 



Pr 



i||UAK+af Hi < 4||AK+af Hi 



k 



<exp(--(r' 



(l-ln/3)) 



< exp 

< exp 



k 



(ln/3-i; 



2 

kin (3 
2~ 



k 

13" 



(2.6) 



By a suitable choice of constants in k = 0(log^ (1/^)) this can be made less than 5. ■ 
Bounding the total of fj, as well as the running time gives the following guarantees about ApproxStr: 

Lemma 2.8 Let A be a n x d matrix, and B an approximation with n rows such that K = B is a 
K-stretch base for A. Then ApproxStr{ A, B, /3, 5) runs in 0{{nnz{A) + nnz{A) + d'^)logij (1/5) + {h + 
d)d^^^ log de~^) time returns upper bounds fj such that: 

1. For each i, with probability at least 1 — 5 we have fj > Tj 

2- \\t\\i < 0{l3Kdlog{l/8)) with probability at least 1 — 5. 

Proof Lemma 2.7 gives 



lUAK^af Hi > \sTn^{a 



/3 
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This and Part 1 of Definition 2.4 means scaling the estimates up by a factor of /3 suffices to give upper 
bounds of the leverage scores. The upper bound on Y^ - fj then can be obtained similar to the proof of 
Lemma 2.7, along with the bound on STTZ-^{mata) given in Part 2 of Lemma 2.6. 
Note that: 

j:f.=f llUAKVlli 

i 

=^||UAK+A^I|2 (2.7) 
k 

Let y be the vector that's the row norms of BK^A. Note that l|y[|2 = X^^Li l|BK^a?"||2 
= 5T7^j^(A) < Kd, where the last inequality follows from Part 2 of Definition 2.4. Part 2 of Lemma 2.6 
can be used to bound ^||Uy|||. Note that when i > 10, t — 1 — In t > t/2. If /c is larger than some absolute 
constant, this gives: 



Pr UllUyll^ >/?n<i>log(l/i)||y||^ 
<o(*exp(-|l5iM)))<i p.8) 

It remains to bound the running time of LeverageUpper. Finding K takes 0{nd^^^) time, while 
inverting it takes an additional time. Computing UA can be done in 0(nnz(A)A;) time, and multiplying 
it on the right by to obtain UAK takes 0{kd'^) time. It remains to evaluate UAK a?" for each row 
i. This can be done by summing nnz(aj) length k vectors, giving a total of nnz(A)/c over all n vectors. 
Therefore the total cost for computing the estimates is 0(/c(nnz(A) + nnz(A) + d'^) + (n + d)d^^'^). ■ 

3 Iterated Row Sampling 

In this section we show an algorithm that computes row samplings by computing upper bounds for the 
leverage scores of blocks of A. It starts from the observation given in Fact 2.3, namely that the stretch 
of row i is the L2 norm of the z*^ row of a matrix, Y. An immediate extension of this is that the 
generalized stretch of a block of rows of A equals to the Forbenius norm of the corresponding rows in Y. 
The Johnson-Lindestrauss lemma as stated in Lemma 2.6 in turn implies that taking a (random) linear 
combination of these rows preserves the total Forbenius norm. In other words, when the original matrix 
is used as the base, the generalized stretch of the smaller set of rows w.r.t. A is comparable to that of 
the original block. Somewhat surprisingly, switching the base to the new matrix only introduces only 
a small error. Therefore, the probabilities required by Lemma 2.1 for computing row samplings can be 
obtained on this smaller matrix. Once again, we need to approximate the product of this smaller matrix 
to obtain these estimates. However, the significantly reduced row count means we can iterate on it by 
computing another projection. In a constant number of iterations the number of rows becomes fewer than 
nnz(A)/d, at which point we can directly compute the answer. We then show that approximate leverage 
scores computed on these smaller projections can be prolongated to earlier, larger matrices while still 
meeting the requirements of Lemma 2.1. This leads to a process that resembles V-cycle algorithms from 
the multigrid literature [BHMOO]. These algorithms create a sequence of coarser versions of the initial 
matrix, and compute answers by going once downward, and once upward through the entire sequence of 
problems. A representation of our algorithm is given in Figure 1 

3.1 Projections and Prolongations 

We first formalize the processes of projecting a matrix to one with fewer rows, and prolonging estimates on 
the projection back to the original matrix. Our key operation is to combine every R rows into k rows, where 
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A(0) 



Afll 



B = B(0) B(l) 



1/ 



f(l) 



\ 



ML - 1) 

B(L-l) 
/ ^ 



ML] 



r{L)\ 



Figure 1: High level structure of algorithm: projections to build a sequence of smaller matrices, then 
prolongations to convert solutions for smaller matrices to solutions for larger ones. 

R and k are parameters that we set to and a constant respectively, . By padding A with additional 
rows of zeros, we may assume that the number of rows is divisible by R. We will use rib = n/R to denote 
the number of blocks, and use the notation to index into the b^^ block. Projection and prolongation 
can be formally defined as: 

1. (R, /i;)-projection: for each block A(^), pick U^;,) to be a /c x i? random Gaussian matrix with entries 
picked independently from A/'(0, 1) and compute A{short)(^i,) = ^(b)-^(b)- Concatenate the blocks 
A{short) (^h) together vertically to form A{short). 

2. (k, i?)-prolongation: given a vector f' with length nf,k, maps it to length n^R vector where each entry 
in block b has value ||t^'^)||i. 

We first show that projections preserve the stretch of blocks w.r.t. A. This can be done by bounding 
the effect of U({,) on the norm of each column of A^^^K^^^^. It follows directly properties of Johnson- 
Lindenstrauss projections described in Lemma 2.6, and a proof is given in Appendix B. 

Lemma 3.1 Assume R > for some constant a and let A{short) be a (R, k) -projection of A. For any 
constant c > there exists a constant k = 0{c/a) such that for any block b, we have with probability at 
least 1 — d~'^ : 

STnA{A{short)^b)) >^STnA{A^b)) 

We next show that we can change the base from A to A{short), and use <ST7^A{s/iort)(-A-(s/iort)(j)) as 
upper bounds for STTZA[A{short)(i,))- As a first step, we need to relate K to the product of A{short), 
K.{short) = A{short)^ A{short) . Since each A{shori)(i,) is formed by merging rows of Aj-^), A{short)^b) — 
U(f,)A({,), A{short)J^^A{short)Qj) can be upper bounded by ||U(6)|||. • A^-)A(b). 

Lemma 3.2 The following holds for each block b: 

A{short)Jf^^A{short)^h) ^| I I If • ^5)^(6) 

A proof is given in Appendix B. Generalized stretches w.r.t. K.{short) and K are evaluated under the 
norms given by K.{short)'^ and K^. Therefore, we need the inverse of this bound, which we obtain using 
the following lemma about pseudoinverses. 

Lemma 3.3 Let C and D be symmetric positive semi-definite matrices and let 11 be the projection space 
onto the range space of C. Then: 

nc^n >z n(c+ z))+n 
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This statement is straightforward when both C and D are full rank, or share the same null space. 
However, as pseudo-inverses do not act on the null space, it is crucial that we're only considering vectors 
of the form a[. This Lemma is also proven in Appendix B. Combining it with bounds on HUi^f,)!!!, allows 
us to bound the distortion caused by switching the base matrix from A to A(short). 

Lemma 3.4 For any constant c, there exists a constant d such that with probability at most 1 — for 
any row i of A{short), a{short)i we have: 

c'lognkR • STllA{short){o.{short)i) > STlZA{short){a{short)i) 

Proof Since each U(6) consists of k x R independent random variables chosen from A^(0, 1), ||U(6)|||, is 
distributed as A^(0, kR). This gives: 

Pr{\\V{b)\\l>l-kR) <exp{-l) (3.9) 

There exists c' such that setting I = c'logn allows us to upper bound this probability by Taking a 

union bound over all n/d blocks gives that with probability at least 1 — d'"^ this holds for all blocks. 

As before, let K = A"^A and K(s/iori) = A{short)^ A{short) . Applying Lemma 3.2 to all blocks and 
summing over all blocks gives: 

K{short) < d log nkR ■ K (3.10) 

Let n be the projection operator onto the range space of K(s/iort), Applying Lemma 3.3 with C = K(s/iort) 
and D = c' log nkR ■ K gives: 

nK(s/iort)+n ^— ^— — HK+n (3.11) 

c log nkR 

Note that SL^short)^ is completely contained within the range space of K.{short). Therefore Ila.{short)f = 
a{short)f and: 

dlognkR ■ ST'JlA{shart)i^{short)i) 
=c' lognkR ■ a{short)i'K{short)^ aj 
=d lognkR ■ a{short)iIl'K{short)~^Ila{short)J 
>a{short)iIl'K^Ila{short)J By Lemma 3.3 

=a{short)iK+a{short)J = STnA{ai{short)i) (3.12) 

Summing over all k rows of A{short)(p-^ completes the proof. I 
Combining Lemmas 3.4 and 3.1 shows that with high probability, scaling up '57~7^A(s/iort) 

(A(s/iort)(b)) 

by 0(ii^ log n) gives upper bounds for a (&(&)) • 

Corollary 3.5 For any parameter R > d", constant c, there exist constants k and d such that for a block 
b, with probability at least 1 — d"'^, 

dR^ log nSTnA{short){Mshort) (J,)) > STTlAia^p)) 

3.2 Iterative Algorithm 

Projecting A to A(short) gives a matrix with fewer rows, and a way to reduce the sizes of our problems. 
A fast algorithm follows by examining the sequence of matrices A = A(0), A(l), . . . A(L) obtained using 
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such projections. Once A(L) has fewer than nnz(A)(i~^ rows, A(L)-^A(L) can be approximated directly. 
This then ahows us to approximate the statistical leverage scores of the rows of A(L) Corollary 3.5 shows 
that statistical leverage scores in A(/), r(/) can upper bound those of A(Z — 1), t{1 — 1). This means we 
can gradually prolongate solutions backwards from A(L) to A(0). We do so by maintaining the invariant 
that B(/) has a small number of rows and is close to A(^). This means stretches of A(/) w.r.t. B(Z) as base 
can be used in place of statistical leverage scores of A(/ — 1). These scores can then be used to compute 
B(/ — 1), keeping the invariant for / — 1. Pseudocode of the algorithm is shown in Algorithm 2. 

Algorithm 2 Row Sampling Via. Projections 
RowCombine(A, R, e, 6, nmin) 

Input: Reduction rate R, n x d matrix A, allowed approximation error e, failure probability 5 = d~'^. 
Output: Sparsifier B that contains 0{dR^\ogd/ e^) scaled rows of A such that (1 - e)A^A < B^B < 
(1 + e)A^A. 

Set L = d^ logjij, and fix constants c', k based on c 

Set e(0) =e, e(l) . . . e(0 = 1/2 

A(0) = A 

for / = 1 . . . L do 

Let A{1) be a (i?, /c)-projection of A(/ — 1) 
end for 
B(L) ^ A(L) 
for i = L . . . 1 do 

f'{l) ^ ^ + 0{R^\ogn) ■ AppROxSTR(A(/),y|B(/),fl,d-^-6) 

Let f{l — 1) be a (fc, i?)-prolongation of f '(/) 

B(/ - 1) ^ Sample(A(/ - l),f(/ - l),e(0) 
end for 
return B^*^^ 



The algorithm obtains sampling probabilities for A(/ — 1) using {k,R) prolongations of estimates of 
stretches of A(/) w.r.t. B(/). We first show these values, f{l — 1), are with high probability upper bounds 
for the statistical leverage scores of A(/ — 1), r(/ — 1). 

Lemma 3.6 Suppose B{1) satisfies: 

lA{lfA{l) r< B{lfB{l) ^ ^A{lfA{l) 

Then with probability at least 1 — 2d^'^^^, we have 

• f{l - 1) > t{1 - 1) 

• I|f(/-l)I|i <0(dii^logn + ^) 

Proof By Lemma 2.5, ^J^B{l) is a 3d-stretch base for A(/ + 1). Since A(/) is a projection of A(/ — 1), 
we can index corresponding blocks in them. Consider any block (6). Applying Lemma 2.8 to all -R < d 
rows of A(^)(b) gives that we have with probability at least 1 — d~'^~^: 

i2||r'(/)(,)||i >5r7^A(/)(A(0(^,)) (3.13) 
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Combining with Corollary 3.5 applied to A(/ — 1) and A(/) gives that with probability at least l — 2d ^ ^ 
we have: 

c'i?3iogn||f'(/)(,)||i > \\t{1 - l)(fe)||i (3.14) 

Also, note that ^ is added to all entries of f (/). Since ||t(/)||i < d, the total number of blocks where ^ 
does not suffice as an upper bound for some row is at most d'^. Applying the above bound to these blocks 
gives that with probability of at least 1 — d^^^"^ we have f(/ — 1) > r(/ — 1). 

It remains to upper bound ||f(/ — l)||i. Lemma 2.8 Part 2 gives with probability at least 1 — d~^~^, 
< 0{d). As each ||T'(/)(f,)||i is prolonged to R rows of f{l — 1), we get: 



if (/ -1)\\^=Y^r(1^ + C- i?3logn||f' (0(6)||i) 
b ^ ' 

<0 (^:^+R^lognY,\\r'{%)\\?j 



<0(dR^logn + ^j (3.15) 

The overall probability of 1 — d^^^^ follows applying union bound once again. ■ 
Combining these with the fact that the number of rows decrease by a factor of 0{R) per iteration 
completes the algorithm. The bound given in Theorem 1.1 is then obtained by setting R to d°'/'^. Applying 
3.6 inductively backwards in I gives the overall bound, and in turn Theorem 1.1. 

Lemma 3.7 /f RowCombine is ran with R = d"/^, 5 = d~'^ where a and c are positive constants, then 
with probability at least 1 — d~^ it returns in 0{nnz{A) + d^~^'^e~'^) time B such that: 

(1 - e)A^A < B^B + e)A^A 

and B has 0{{d^^°^ logn + ^) logde"^) rows, each being a scaled copy of some row of A, 

Proof We first show correctness via. induction backwards on /. Specifically we show that with probability 
at least 1 - 3(L - ^)(i"'="^ B(0 has 0{{dR^ log n + ^) log de{l)~^) rows and satisfies (1 - e(/))A(O^A(/)^ < 
B(I)'^B(/) ^ (1 + e(/))A(/)^A(/)^. Note that since R = d"' , L = log^d^ is a constant. As the number of 
rows decreases by a factor of 0{R) per projection, we get that A(L) has 0{n/d^) rows. Therefore the base 
case where I = L follows from B(L) = A(L). 

For the inductive step, we assume that the inductive hypothesis holds for / > 1 and try to show it for 
I — 1. As e(Z) was set to 1/2, we have: 

\A{ifA{i) r< mfm ^ Inifm 

This allows us to invoke Lemma 3.6, which combined with Lemma 2.1 gives that with probability 1 — d^^^^, 
the inductive hypothesis also holds for l — l. The failure probability follows from union bounding this with 
the failure probabilities from the hypothesis and Lemma 3.6. 

We now bound the total running time, starting with the projections. Since '^(b) has k = 0(1) rows, 
the total cost of computing U(;,) A(/)(b) is 0(nnz(A(Z))). Furthermore, we also get nnz(A(Z + 1)) = 
0(nnz(A(/))), giving a total cost of 0(nnz(A)) over all L projections. Since B(/) has 0{{dR'^logn + 
^)logde{l + 1)~^) rows. Lemma 2.8 gives that each cah to ApproxStr takes 0((i'^ii'^ log n log de~^ + 
nd^^'^ log de^"^) time. By standard assumptions that logde"^ < d^~^ and n < nnz(A), this is at most 
0(d'^i?^ log 71 log (ie~^ + nnz(A)). Substituting R = d"/^ into these bounds gives the result. ■ 
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4 Introduction and Removal of Artificial Bases 

This section describes a different algorithm that computes row samplings of A using artificial bases. These 
bases are formed by adding d additional rows to A equaling to 7I. When 7 very large, this matrix can 
be shown to be almost the same as the identity matrix. On the other hand, when 7 is very small, there 
is very little difference between this matrix and A. As a result, the key is to show that 7 can be reduced 
geometrically over the iterations. Given an initial value 7(0), and a parameter R indicating the reduction 
rate, we consider a geometrically decreasing sequence of 7(/)s defined by 7(/ + 1) = R~^^'^'y(l). Then A(Z) 
can be written as: 



A(/) 



A 

7(01 



And we will use K(Z) to denote its product: 

K(/) = A(/f A(Z) = K + 7(Z)2l 

We will use B(/) to denote approximations of A(/) obtained from row samplings, and K(l) to be the 
product of B(Z), B(Z)^B(Z). A geometrically decreasing sequence of values of 7 are of special interest to us. 
We first show that if 7(0)^ is more than the maximum eigenvalue of A^A, a copy of the identity matrix 
suffices as B(0). 

Lemma 4.1 Let Amax be the maximum eigenvalue of A^A. Ifj(0)^ > Amax; then B(l) = 7(0)/ satisfies: 

k(0) ^ K{0) ^ 2k{0) 

Proof The LHS follows from K(/) = K + -f{Ofl h 7(0)^1 = B(0)'^B(0). The RHS can be obtained by 
K(0) y AmaxI >1 7(0)^1 h K(0) and adding -f{Ofl to both sides. ■ 
Also, once 7(0^ is less than the minimum eigenvalue of A'^A, B(Z) can be directly used as a low stretch 
base for A(Z). 



Lemma 4.2 Let Amin be the minimum eigenvalue of A. If ^(Xf' < Xmm, then y^B{l) is a 2d-stretch 
base for A. 



Let a spectral decomposition of K be K = Ajvv^, where Ai < A2 < . . . < A,^. Assume Afe is the first 



Proof For this proof it is easier to work with the product of ^/^B(Z), since (y ^B(/))-^(-^/ ^B(Z)) = ^K(Z). 

Let a spectral decomposite 
non-zero eigenvalue, then: 

K(/) =K + 7(0'I = ^(A, + 7(/)')v,v] 



d 

3 

j=k 



Consider the eigenvector Vj for some j < k, Xj = implies ||Avi||| = v^Kv = 0, Therefore for any 
row i of A, we have a^Vj = 0. This means for both ajK+a?' and ajK(Z)+a?^ it suffices to consider the 
terms from k to d. For each of these terms, 7 < Xj gives ^(Aj + 7) < Xj and therefore: 

a; - (l/2)(A, +7) ^"^-^^^ 
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As each term from k to d is non-negative, summing over them meets Part 1 of Definition 2.4. 
For the upper bound required by Part 2, note that Xj < Aj ■ + 7 gives: 

1 2 

<- 



(l/2)(A,+7) -A,- 
For a row i of A, a^, summing over each eigenvalue gives: 

. Summing once again over all rows leads to STIZ rf (A) < 25T7^a(A) < 2d. ■ 

V 2^(') 

The algorithm is based on the fact that K(/ + 1) can be used as 0(di?)-stretch base for K(/). By the 
requirements of Lemma 2.5, it suffices to show that K(/) and K(/ + 1) are spectrally close. 

Lemma 4.3 

K{1 + K{1) < RK{1 + 1) 

Proof The LHS then follows from I being PSD, while for the RHS we have: 

RK{l + l) =i?(K + 7(/ + l)2l) 
=R(K + R-^^{lfl) 
>=K + 7(0I 

=m (4.17) 

■ 

Note that if ^K{1) ^ K(/) ^ |K(/), Lemma 4.3 implies that 

^K(/ + 1) ^^K(0 ^ K(/ + 1) (4.18) 

Therefore by Lemma 2.5, ^J^^^{l) is a 3(ii?-stretch base for K(/ + 1). If we're given bounds on non-zero 

eigenvalues Xupper and Xiower, we can start with 7(0) = Xupper and continue until ^{l) < Xiower- This leads 
to the algorithm whose pseudocode is given in Algorithm 3 



Theorem 4.4 Let c be any constant corresponding to a failure probability of d ^. Given a n x d matrix 

^upper 



A, upper and lower bounds on the non-zero eigenvalues of A Xupper o-nd Xiower ! (^s well as parameters e. 



RowSample runs in 0(((nnz(A) -|- d^) logj^ n + d^li^ log^ d) logji{Xupper/ Xiower) + d'^R'^ log^ de"^) time, 
and returns a matrix B that contains 0{dR^ log de~'^) rescaled rows of A and K = B^B. With probability 
at least 1 — d~^ we have: 

{l-e)K^k^ il + e)K 

Proof We may assume that log R{Xupper / Xiower) < l/lOd since otherwise directly computing A"^A gives 
the required runtime bound. 

We first show that K(/) are good approximations of K(/) with high probability. We prove by induction 
on / that with probability 1 — 0{ld~'^~^) we have: 

1k(0 ^K(0 ^ ^K(0 (4.19) 
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Algorithm 3 Row Sampling via. Introduction and Removal of Identity 

ROWSAMPLE(i?, Xlower, Kpper, A, €, 5) 

Input: Parameter i? > 4 indicating the speed of reduction. Matrix A, upper and lower bounds on eigen- 
values, Xupper > Amax(A) and Xlower < Amin(A). Error parameter e and failure probability 5. 
Output: Matrix B containing 0{dR'^ logd) scaled rows of A such that (1 - e) A^A ^ B^B ^ (1 + e) A^A 



1 

2 
3 

4 
5 

6 

7 
8 
9 

10 
11 



Let 7(0) = Xu pper 

Let B(0) = Xupperl, K(0) = Xupperl, I = 

while 7(/) > Xlower do 

7(0 ^ R~^'Hi - 1) 

f (/) ^ ApproxStr(A(/), yj'B(/ - 1), R, d-^n-^) 

B(/) ^ Sample(A(/), f (/), 1/2, d-^5) 
end while 

f ^ ApproxStr(A, ^J^B{l), R, d'^'n'^) 

B ^ Sample(A, f, e, 
return B 



and K(/) contains 0{dR?\o^ d) scaled rows of A(Z). The base case of Z = follows from Lemma 4.1. For 
the inductive case, combining Lemma 4.3 and the inductive hypothesis gives: 

^K(/) ^ Ak(/ - 1) ^ K(/) (4.20) 

Therefore Lemma 2.5 gives that ^B(/ — 1) is an 3(i-R-stretch base for A(/). Lemma 2.8 then guarantees 
that for each row i of A(Z), with probability at least 1 — d'^~^{n + d)~^, f{l) > t{1). Applying union bound 
over all n + d rows of A(l) gives that with probability at least 1 — d~'^~^ this holds for all rows. Also, with 
probability at least 1 - d'""-^, \\f{l)\\i < 0{dR^ log(d-^-i)) = 0{dR^ \ogd). If both of these bounds hold, 
by Lemma 2.1, we have with probability at least 1 — 0{d~'^~'^), ^K(Z) ^ K ^ |K(/), and B(/) consists of 
0{dR^ log^ d) rows of A(/). Combining this with the error probability in the induction hypothesis shows 
that it holds for / as well. 

For final output, combining < Amin with Lemma 4.2 gives that y^B(/) is a 2(i-stretch base for A. 

Applying Lemmas 2.8 and Lemma 2.1 once more to f gives the bound on K. 

The last call to Sample produces a B with O {dR"^ log^ d/e"^) scaled rows of A, which means it costs 
0{d^R^log^d/e^) to evaluate K. It remains to bound the cost of approximating stretches and sampling 
from the O (log ji{Xupper / Xiower) iterations. As B(/) has 0{dR^ log^ d) rows. Lemma 2.8 gives a total of 
0((nnz(74) + nnz(B) + d^) log^ n + d'^R'^ log^ d) time. Combining this terms with the iteration count gives 
the overall running time. ■ 

In order to obtain a running time of 0(nnz(A) + for some constant a > 0, we can set R to d"^/^. 

Theorem 4.4 gives that RowSample runs in 0(((nnz(A) + d^) log^ n + log^ d) log^iXupper / Xiower) + 
(]i^+a jQg2 i[YD.e, and with probability at least 1 — 6 returns B with 0{d^~^°' log^ de"^) rescaled rows of 

A. When n and Xupper / Xiower are both poly((i), the running time simplifies to 0(nnz(A) + d^"'"°e~^). 

The pseudocode as given in Algorithm 3 needs values for Xupper and Xiower given as input. This 
requirement can be relaxed in several ways. Taking the Forbenius norm of A, ||A||2 gives a factor d 
approximation of Amax- Furthermore, if we know the rank of A, we can check whether 7(/) < Amin by 
computing the eigenvalues of K(/) and comparing against the smallest one. 
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5 Remarks 

Aside from giving the state-of-art running times for approximating A^A, each row in the output of our 
algorithms is a scaled copy of some row of A. This means any specialized structure for rows of A are 
preserved in B. The algorithm described in Section 4 has the further advantage that all the intermediate 
matrices also consists of only rows of A and a minor of the identity matrix. In practice, this means any 
uniform sparsity in A would also exist in these intermediate problems. In situations where all columns of 
A is sparse, iterative methods can solve the d x d linear systems in times faster than d'^ . This speedup is 
difficult to quantify without the introduction of additional parameters, and there are few general bounds 
for iterative methods with sparse systems. As a result, we believe this additional property is of more 
practical value and should be examined experimentally. 

We believe that removing the factor of d", as well as extending these approaches computing low-rank 
approximations are natural directions for future work. 
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A Proofs Prom Section 2 
Proof of Fact 2.2: 



i=l 



i=l 

=tr 



a, a,; 



=tr [K+K] 



(1.21) 



Proof of Fact 2.3: Note that both the generahzed definition of leverage scores and the Forbenius 
norm works on each row independently. Therefore it suffices to prove this when A' has a single row, aka. 
A' = a. In this case the cyclic property of trace gives: 



5r7^K(a) =aK+a^ 

=aK+V2K+i/2aT 

= ||K+i/2a'r||2 



(1.22) 



Proof of Lemma 2.5: The condition given implies that the null spaces of K and K are identical, 
giving: 



(1.23) 
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Applying this to the vector x gives: 

xK+x^ <xK"*'x'^ < KxK'^'x'^ (1.24) 

The first inequahty gives Part 1, while summing the second inequahty over all rows gives Part 2. ■ 
B Proofs About Projections and Prolongations 

Proof of Lemma 3.1: Let K = A^A, by Fact 2.3, we have that 

5r7eA(A(5)) = ||A(6)K+i/2||2, 

and 

STnA{A{short)(^b)) = l|A(s/iort)(f,)K+^/2||2^ 
Furthermore, since A(s/iort)(b) = U(-ft)A(5), we have: 

STnA{A{short)(^h)) =l|A(s/iort)(b)K+^/^|||. 

=||U(b)A(,)K+i/2|[2, (2.25) 

Let Y denote A(6)K+^/2^ ^^i^^ 

we have: 

Pr (^STnAiA{short)^b)) < ^STTlAiA^b))^ 

=Prr||U(5)Y|||<^||Y|||') (2.26) 



Consider each column of Y, y. j. By Lemma 2.6, we have that: 



Pr(^l|U(5)y,,||i<|lly:,,Ili 
<exp(^-^(/3-(l + ln/3))) (2.27) 

Setting /3 = R, this value can be bounded by 0(A;) exp(-|(/3-i - 1 - In /3) < 0{k)R-^/^. As R > (P, 
setting k = 4(c + l)/a gives that this probability is at most d~'^~^ . Applying union bound over all d 
columns gives that with probability 1 — we have that for all columns, ||U(f,)y. > ^Hyijlll- The 
overall bound then follows from fact that ||Y|||, = ||y jlli- ^ 

Proof of Lemma 3.2: By definition of spectral norm, we have ||U^^U(ft)|| < ||U(6)|||., and in turn: 

l|U(b)y||i <||U(fe)|||.-||y||i (2.28) 

for any vector y. 
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Let X be any vector in 3?"^. Substituting in A(s/iort)({,) = U(-fe)A(-^) gives: 



X A(s/iort)^^-j A(s/iort)(f,)X 
||A(s/iort)(b)x||| 



<I|U, 
=I|U 



(6) I If 

(fe)llF 



•x^Aj)A(6) 



By Equation 2.28 



X 



(2.29) 



Proof of Lemma 3.3: Consider an orthonormal basis for the range space of C, vi . . . ^rank{C)- Since 
C + D >z C, this basis can be extended to an orthonormal basis to the range space of C + D by adding 
Vranfc(c)+i • • ■^rank{c+'D)- It suffices to prove the claim under this basis system. Here C and D can be 
rewritten as: 



Cii 









D 



D 
D 



11 



12 



Dl2 

D22 



(2.30) 



Where Cn and D22 are strictly positive definite. Furthermore, since D is positive semi-definite we have 
that Dii — Di2D^2^Di2 is also positive semi-definite. For any vector x, IIx gives a vector that's non-zero 
only in the first rank{C) entries. Let this part be xi. Then evaluating (C + D)+nx = [yi;y2] becomes 
solving the following system: 



(Cii + Dii)yi + Di2y2 =xi 
Df2yi + D22y2 =0 

The second equation gives y2 = — D22^^D^2yi- Substituting it into the first one gives: 

(Cn + Dn - Di2D22^Df2)yi = xi 



(2.31) 
(2.32) 



(2.33) 



Note that this is the same as taking the partial Cholesky factorization onto the range space of H. Combining 
things gives: 



x'^n(C + D)+nx =xf (Cn + Dn - Di2D22^D2 - 1, 12'^)~^xi 



(2.34) 



Since Dn — Di2D22^D{2 is positive definite, we have Cn ^ C + Dn — Di2D22^D{2 and therefore: 



-ii 



x^n(C + D)+nx =xf (Cn + Dn 
^x^ "-^n xi 

=x^nc+nx 



Di2D2-2^Df2)xf 



(2.35) 
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